/**
 * @file   test.cpp
 * @author Kimmo Luoma <kimmo@PHOTON3>
 * @date   Tue Oct  9 12:43:05 2012
 * 
 * @brief  Simple test for coin, state and shift class
 * 
 * 
 */

#include "coin.hpp"
#include "shift.hpp"
#include "state.hpp"
#include<boost/math/constants/constants.hpp>

int main(){
  
  //We use ring example
  int N=9;
  int d=2;
  int T=100;
  
  matrix<std::complex<double> > psi(d,N,0);
  psi(0,0)=1/std::sqrt(2);
  psi(1,0)=std::complex<double>(0,1/std::sqrt(2));
  
  //DFT coin
  std::complex<double> I(0,1);
  const double pi=boost::math::constants::pi<double>();
  matrix<std::complex<double> > H(d,d);
  for(int k=0;k < d; k++){
    for(int l=0; l < d; l++){
      H(k,l)=std::exp(2*pi*I*k*l/d);
    }
  }
  H=H*(1/std::sqrt(d));
  
  /**
   * Test all constructors for coin op
   * 
   */
  QuantumWalk::coin<std::complex<double > > C0;
  QuantumWalk::coin<std::complex<double > > CH(H,N,d);
  C0.set_op(H);
  QuantumWalk::coin<std::complex<double > > C1(CH);
  QuantumWalk::coin<std::complex<double > > C2=CH;

  std::cout<<"C0.get_op() "<<C0.get_op()<<std::endl;
  std::cout<<"C1.get_op() "<<C1.get_op()<<std::endl;
  std::cout<<"C2.get_op() "<<C2.get_op()<<std::endl;
  std::cout<<"CH(0,0) "<<CH(0,0)<<std::endl;
  

  /**
   * Teststate. Example is ring
   * 
   */
  /// Default
  QuantumWalk::state<std::complex<double> > s;
  QuantumWalk::state<std::complex<double> > s1(psi);
  QuantumWalk::state<std::complex<double> > s2(s1);
  s=s1;
  std::cout<<"Site probabilities s "<<s.prob_sites()<<std::endl;
  std::cout<<"Internal probabilities s"<<s.prob_internal()<<std::endl;
  std::cout<<"Site probabilities s1 "<<s1.prob_sites()<<std::endl;
  std::cout<<"Internal probabilities s1"<<s1.prob_internal()<<std::endl;
  std::cout<<"Site probabilities s2 "<<s2.prob_sites()<<std::endl;
  std::cout<<"Internal probabilities s2"<<s2.prob_internal()<<std::endl;
  
  /**
   * Test shift operator
   * 
   */
  matrix<std::complex<double> > p1(N,N,0),p2(N,N,0);
  for(int k=0;k<N;k++){
    p1((k+1) % N,k)=1;
    p2(k,(k+1) % N)=1;
  }
  QuantumWalk::shift<std::complex<double> > S;
  S.add_op(0,p1);
  S.add_op(1,p2);
  
  
  for(int k=0;k<T;k++){
    s*=CH;
    s*=S;
  }
  std::cout<<"State after "<<T<<" steps:"<<s.prob_sites()<<std::endl;
  std::cout<<s.prob_internal()<<std::endl;
  return 0;
}
